 %% Load data
load /Users/rogerfu/Dropbox/QDM2_data/Speleothems/Onca/Limb_data_RF/compiled_images/BC_Limb_Bz_cal.mat
load /Users/rogerfu/Dropbox/QDM2_data/Speleothems/Onca/Limb_data_RF/compiled_images/AB_Limb_Bz_cal.mat
load /Users/rogerfu/Dropbox/QDM2_data/Speleothems/Onca/Limb_data_RF/compiled_images/CD_Limb_Bz_cal.mat
load /Users/rogerfu/Dropbox/QDM2_data/Speleothems/Onca/Limb_data_RF/compiled_images/DE_Limb_Bz_cal.mat
load('/Users/rogerfu/Dropbox/QDM2_data/Speleothems/Onca/Paleoclimate_class_RF/yearly_compiled_Bz.mat')
load /Users/rogerfu/Dropbox/QDM2_data/Speleothems/Onca/Limb_data_RF/compiled_images/annual_precip.mat
load /Users/rogerfu/Dropbox/QDM2_data/Speleothems/Onca/O_18_annual_onca2_1913.mat
load /Users/rogerfu/Dropbox/QDM2_data/Speleothems/Onca/Carbon_13_annual_1913.mat
%% Load data - KIM
close all
clc
clear all
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/Limb_data_RF/compiled_images/BC_Limb_Bz_cal.mat');
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/Limb_data_RF/compiled_images/AB_Limb_Bz_cal.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/Limb_data_RF/compiled_images/CD_Limb_Bz_cal.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/Limb_data_RF/compiled_images/DE_Limb_Bz_cal.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/Paleoclimate_class_RF/yearly_compiled_Bz.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/Limb_data_RF/compiled_images/annual_precip.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/O_18_annual_onca2_1913.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/Carbon_13_annual_1913.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/rain_total.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/binned_Bz.mat')
load ('/Users/kimberlyhess/Harvard Paleomag Dropbox/Harvard Paleomag/QDM2_data/Speleothems/ONCA/Onca/rain_total.mat')
%%
A=5

%% plotting unsmoothed data and trimming anomalies near speleothem surface
%remove some extreme values early on in the time series
AB_Limb_Bz=AB_Limb_Bz(AB_Limb_Bz(:,5)>1994,:);
BC_Limb_Bz=BC_Limb_Bz(BC_Limb_Bz(:,5)<2016,:);

%CD_Limb_Bz=CD_Limb_Bz(CD_Limb_Bz(:,5)<1990,:);
CD_Limb_Bz=CD_Limb_Bz(CD_Limb_Bz(:,5)>1951,:);
DE_Limb_Bz=DE_Limb_Bz(DE_Limb_Bz(:,5)>1913,:);
yearly_compiled_Bz=yearly_compiled_Bz(yearly_compiled_Bz(:,1)<2016,:);


plot_no_smoothing=true;

rainscale=max(-yearly_compiled_Bz(:,2))/...
    max(annual_precip(:,2));


if plot_no_smoothing
    figure
    plot (AB_Limb_Bz(:, 5), AB_Limb_Bz(:,2),'r')
    hold on
    plot (BC_Limb_Bz(:, 5), BC_Limb_Bz(:,2),'b')
    hold on
    plot (CD_Limb_Bz(:, 5), CD_Limb_Bz(:,2),'g')
    hold on
    plot (DE_Limb_Bz(:, 5), DE_Limb_Bz(:,2),'k')
    hold on
%     plot (yearly_compiled_Bz(:, 1), -yearly_compiled_Bz(:,2),'k--')
    hold on
   % plot (annual_precip(:, 1), rainscale*annual_precip(:,2),'b--')
end


%% making time series arrays 
%these are just a list of time coordinates for each profile

step=1; %in years


% one for each limb profile
time_seriesAB=double(floor(min(AB_Limb_Bz(:,5))):step:ceil(max(AB_Limb_Bz(:,5))));
time_seriesBC=floor(min(BC_Limb_Bz(:,5))):step:ceil(max(BC_Limb_Bz(:,5)));
time_seriesCD=floor(min(CD_Limb_Bz(:,5))):step:ceil(max(CD_Limb_Bz(:,5)));
time_seriesDE=floor(min(DE_Limb_Bz(:,5))):step:ceil(max(DE_Limb_Bz(:,5)));
% for the center transect
time_seriesCENTER=floor(min(yearly_compiled_Bz(:,1))):step:ceil(max(yearly_compiled_Bz(:,1)));
% for the precip data
time_seriesRAIN=floor(min(annual_precip(:,1))):step:ceil(max(annual_precip(:,1)));
% for the O18 data
time_seriesOXYGEN=floor(min(Oxygen_18_annual(:,1))):step:ceil(max(Oxygen_18_annual(:,1))); 
% for the C13 data
time_seriesCARBON=floor(min(Carbon_13_annual(:,1))):step:ceil(max(Carbon_13_annual(:,1))); 


%% applying spline smoothing and plotting 
% does not need to be run to analyze the combined limb profile

plot_smoothed = true;

spline_smoothing_p=.12; %this is the important spline smoothing parameter. smaller = more smoothing

mag_interpBC=csaps(BC_Limb_Bz(:,5),BC_Limb_Bz(:,2),spline_smoothing_p,time_seriesBC);

test=double(AB_Limb_Bz);
mag_interpAB=csaps(test(:,5),test(:,2),spline_smoothing_p,time_seriesAB);

mag_interpCD=csaps(CD_Limb_Bz(:,5),CD_Limb_Bz(:,2),spline_smoothing_p,time_seriesCD);

mag_interpDE=csaps(DE_Limb_Bz(:,5),DE_Limb_Bz(:,2),spline_smoothing_p,time_seriesDE);

mag_interpCENTER=csaps(yearly_compiled_Bz(:,1),yearly_compiled_Bz(:,2),spline_smoothing_p,time_seriesCENTER);

rain_interpRAIN2=csaps(annual_precip(:,1),annual_precip(:,2),spline_smoothing_p,time_seriesRAIN);

Oxygen_interp=csaps(Oxygen_18_annual(:,1),Oxygen_18_annual(:,2),spline_smoothing_p,time_seriesOXYGEN);

Carbon_interp=csaps(Carbon_13_annual(:,1),Carbon_13_annual(:,2),spline_smoothing_p,time_seriesCARBON);

if plot_smoothed
    figure
    plot(time_seriesBC, mag_interpBC)
    hold on
    plot(time_seriesAB, mag_interpAB,'r')
    hold on
    plot(time_seriesCD, mag_interpCD,'g')
    hold on
    plot(time_seriesDE, mag_interpDE,'m')
    hold on
%     plot(time_seriesCENTER, -1*mag_interpCENTER,'k--')
%     hold on
%     plot(time_seriesRAIN, rainscale*rain_interpRAIN2,'b--')
%     hold on
%     plot(time_seriesOXYGEN, Oxygen_interp,'r--')
%     hold on 
%     plot(time_seriesCARBON, Carbon_interp,'g--')
end

%% Averaging the limbs to create a single limb profile

combined_time_series=(min([time_seriesAB time_seriesBC time_seriesCD time_seriesDE]):1:...
    max([time_seriesAB time_seriesBC time_seriesCD time_seriesDE]))';
alllimbs={AB_Limb_Bz, BC_Limb_Bz, CD_Limb_Bz, DE_Limb_Bz};

%average value if the limb series has the year
averagelimbs=zeros(size(combined_time_series,1),2);
for i=1:size(combined_time_series,1) %i runs across combined years
    temp=[];
    for j=1:size(alllimbs,2) %j runs across the limbs
        tempcell=alllimbs{j};
        temp2=[];
        for k=1:size(tempcell,1) %k runs across limb years
            if round(tempcell(k,5)) == combined_time_series(i)
                temp2=[temp2 tempcell(k,2)];
            end
        end
        if isnan(mean(temp2))
        else
            temp=[temp mean(temp2)];
        end
    end
    averagelimbs(i,:)=[combined_time_series(i) mean(temp)];
end

figure
plot (averagelimbs(:, 1), averagelimbs(:,2))
%writematrix(averagelimbs(2:end-1,1:2),'/Users/rogerfu/Dropbox/my_classes/EPS_102_2022/Lectures/11_smoothing/mag_data.csv');
%writematrix([],'/Users/rogerfu/Dropbox/my_classes/EPS_102_2022/Lectures/11_smoothing/rain.csv');

%% averaged limbs with smoothing
%be aware these two parameters are also used earlier, but are defined independantly
step=1; %in years
spline_smoothing_p=.12; %lower is smoother

time_seriesALLLIMBS=floor(min(averagelimbs(:,1))):step:ceil(max(averagelimbs(:,1)));
mag_interpALLLIMBS=csaps(averagelimbs(:,1),averagelimbs(:,2),spline_smoothing_p,time_seriesALLLIMBS);

% time_seriesCENTER=floor(min(yearly_compiled_Bz(:,1))):step:ceil(max(yearly_compiled_Bz(:,1)));
% mag_interpCENTER=csaps(yearly_compiled_Bz(:,1),-yearly_compiled_Bz(:,2),spline_smoothing_p,time_seriesCENTER);
% 
% time_seriesRAIN=floor(min(annual_precip(:,1))):step:ceil(max(annual_precip(:,1)));
% rain_interpRAIN2=csaps(annual_precip(:,1),rainscale*annual_precip(:,2),spline_smoothing_p,time_seriesRAIN);


figure
yyaxis left
plot(time_seriesALLLIMBS, mag_interpALLLIMBS,'r')
hold on
% plot(time_seriesCENTER, mag_interpCENTER,'g')
% hold on
% yyaxis right
% plot(time_seriesRAIN, rain_interpRAIN2,'b--')


%% statistics limb Bz vs rain

figure;

subplot(2, 2, 1); hold on;
title('All data', 'FontSize', 30)
yyaxis left
plot(time_seriesRAIN, rain_interpRAIN2)
ylabel('Yearly rainfall (mm/year)')
yyaxis right
plot(time_seriesALLLIMBS, mag_interpALLLIMBS)
ylabel('Limb Bz')
plot(time_seriesCENTER, mag_interpCENTER,'k')

%make sure time series for mag and rain are of the same length
firstyear=min([time_seriesRAIN time_seriesALLLIMBS time_seriesCENTER]);
lastyear=max([time_seriesRAIN time_seriesALLLIMBS time_seriesCENTER]);
allcombineddata=zeros(lastyear-firstyear+1,4);
allcombineddata(:,1)=firstyear:step:lastyear;
allcombineddata(min(time_seriesRAIN)-firstyear+1:min(time_seriesRAIN)-firstyear+1+size(time_seriesRAIN,2)-1,2)=...
    rain_interpRAIN2/rainscale;
allcombineddata(min(time_seriesALLLIMBS)-firstyear+1:min(time_seriesALLLIMBS)-firstyear+1+size(time_seriesALLLIMBS,2)-1,3)=...
    (mag_interpALLLIMBS);
allcombineddata(min(time_seriesCENTER)-firstyear+1:min(time_seriesCENTER)-firstyear+1+size(time_seriesCENTER,2)-1,4)=...
    (mag_interpCENTER);
overlappyingarray=allcombineddata(:,4).*allcombineddata(:,2).*allcombineddata(:,3);
%column 1 is time; column 2 is rain; column 3 is limb Bz; column 4 is
%center Bz.  All are smoothed. 
alldataoverlapping=allcombineddata(overlappyingarray~=0,:);
%combined with d18O:
alldataoverlapping=[alldataoverlapping Oxygen_interp(1:size(alldataoverlapping,1))'];

columntoregress=4; % 3 = limb vs. rain; 4 = center vs. rain; 5 = 18O vs. rain
yearrange = 78;

subplot(2, 2, 2); hold on;
title(strcat('Smoothed (spline p=', string(spline_smoothing_p), ')'), 'FontSize', 30)
yyaxis left
plot(alldataoverlapping(1:yearrange,1), alldataoverlapping(1:yearrange,2))
ylabel('Yearly rainfall (mm/year)')
yyaxis right
plot(alldataoverlapping(1:yearrange,1), alldataoverlapping(1:yearrange,columntoregress))
ylabel('Bz')

subplot(2, 2, 3); hold on;
title('Regression', 'FontSize', 30)
X = [ones(length(alldataoverlapping(1:yearrange,2)), 1) alldataoverlapping(1:yearrange,2)];
b = X\ alldataoverlapping(1:yearrange,columntoregress);
mag_calc = X*b;
scatter(alldataoverlapping(1:yearrange,2),  alldataoverlapping(1:yearrange,columntoregress))
plot(alldataoverlapping(1:yearrange,2), mag_calc)
xlabel('Yearly rainfall (mm/year)')
ylabel('B_z')
rsq = 1 - sum(( alldataoverlapping(1:yearrange,columntoregress) - mag_calc).^2)/sum(( alldataoverlapping(1:yearrange,columntoregress) -...
    mean( alldataoverlapping(1:yearrange,columntoregress))).^2);
legend(['b = ', num2str(b(2)), ', R^2 = ', num2str(rsq)])

subplot(2, 2, 4); hold on;
title('B_z Residuals', 'FontSize', 30)

plot(alldataoverlapping(1:yearrange,1), alldataoverlapping(1:yearrange,columntoregress)-mag_calc, 'o--')

lm = fitlm(alldataoverlapping(1:yearrange,2), alldataoverlapping(1:yearrange,columntoregress));
lm.Coefficients
lm.Rsquared.Ordinary



%% statistics center vs limb Bz

figure;

subplot(2, 1, 1); hold on;
title('Limb vs. Center', 'FontSize', 30)
yyaxis left
plot(time_seriesALLLIMBS, mag_interpALLLIMBS)
ylabel('Bz')
plot(time_seriesCENTER, mag_interpCENTER,'k')

%make sure time series for mag and rain are of the same length
firstyear=min([time_seriesRAIN time_seriesALLLIMBS time_seriesCENTER]);
lastyear=max([time_seriesRAIN time_seriesALLLIMBS time_seriesCENTER]);
allcombineddata=zeros(lastyear-firstyear+1,4);
allcombineddata(:,1)=firstyear:step:lastyear;
allcombineddata(min(time_seriesRAIN)-firstyear+1:min(time_seriesRAIN)-firstyear+1+size(time_seriesRAIN,2)-1,2)=...
    rain_interpRAIN2;
allcombineddata(min(time_seriesALLLIMBS)-firstyear+1:min(time_seriesALLLIMBS)-firstyear+1+size(time_seriesALLLIMBS,2)-1,3)=...
    (mag_interpALLLIMBS);
allcombineddata(min(time_seriesCENTER)-firstyear+1:min(time_seriesCENTER)-firstyear+1+size(time_seriesCENTER,2)-1,4)=...
    (mag_interpCENTER);
overlappyingarray=allcombineddata(:,4).*allcombineddata(:,2).*allcombineddata(:,3);
%column 1 is time; column 2 is rain; column 3 is limb Bz; column 4 is
%center Bz.  All are smoothed. 
alldataoverlapping=allcombineddata(overlappyingarray~=0,:);

subplot(2, 1, 2); hold on;
title('Regression', 'FontSize', 30)
X = [ones(length(alldataoverlapping(:,3)), 1) alldataoverlapping(:,3)];
b = X\ alldataoverlapping(:,4);
mag_calc = X*b;
scatter(alldataoverlapping(:,3),  alldataoverlapping(:,4))
plot(alldataoverlapping(:,3), mag_calc)
xlabel('Yearly rainfall (mm/year)')
ylabel('B_z')
rsq = 1 - sum(( alldataoverlapping(:,4) - mag_calc).^2)/sum(( alldataoverlapping(:,4) -...
    mean( alldataoverlapping(:,4))).^2);
legend(['b = ', num2str(b(2)), ', R^2 = ', num2str(rsq)])

lm = fitlm(alldataoverlapping(:,3), alldataoverlapping(:,4));
lm.Coefficients

%%
%%

figure (2)
yyaxis left
plot(time_seriesALLLIMBS, mag_interpALLLIMBS)
hold on
plot(time_seriesCENTER, mag_interpCENTER,'k')
hold on
yyaxis right
plot(time_seriesRAIN, rain_interpRAIN2,'b--')

legend

%%

%plot O18 C rain and bz vertically for publication - 4 year
%close all
figure (1)

hold on

%1
subplot(1, 4,1)
y1=(time_seriesALLLIMBS);
x1= (mag_interpALLLIMBS);
plot(x1, y1, 'Color', '#D95319', 'Linewidth', 2)
hold on
y2= (time_seriesCENTER);
x2= (mag_interpCENTER);
set(gca, 'FontSize',15)
plot(x2, y2, 'Color', [0.2 .5 0.5], 'Linewidth', 2)

ylim([1913 2016])
title('Magnetization')
xlabel('Bz')
ylabel('years')

%2
subplot(1,4,2); 
y=(time_seriesRAIN);
x3 = (rain_interpRAIN2/rainscale);
plot(x3,y, 'Color', '#0072BD', 'Linewidth', 2)
set(gca, 'FontSize',15)
title('rain')
ylim([1913 2016])
xlabel('mm rain')
% ylabel('years')

%3
subplot(1,4,3)
y=time_seriesOXYGEN;
x4=Oxygen_interp;
plot(x4, y, 'Color', '#7E2F8E', 'Linewidth', 2)
set(gca, 'FontSize',15)
ylim([1913 2016])
title('{\delta} Oxygen 18')
xlabel('O18')
% ylabel('years')

%4
subplot(1,4,4); 
y=time_seriesCARBON;
x5 = Carbon_interp;
plot(x5,y, 'Color', '#A2142F', 'Linewidth', 2)
set(gca, 'FontSize',15)
title('{\delta} Carbon 13')
ylim([1913 2016])
xlabel('C13')
% ylabel('years')

%% Regression Figures
%first smooth central col BZ data


%smooth rain data

step=1

plot_smoothed = true;

spline_smoothing_p=.12; 
% for the center transect
time_seriesCENTER_Bz=floor(min(yearly_compiled_Bz(:,1))):step:ceil(max(yearly_compiled_Bz(:,1)));

Bz_interp_Cent_final=csaps(yearly_compiled_Bz(:,1),yearly_compiled_Bz(:,2),spline_smoothing_p,time_seriesCENTER_Bz);

Rain_interp_Cent_final=csaps(rain_total(:,1),rain_total(:,2),spline_smoothing_p,time_seriesCENTER_Bz);

%% plotting raw, annual rain to annual bz - central col

figure
x1=binned_rain_total(:,1);
y1= binned_rain_total(:,2);
y2= abs(binned_Bz (:,2));

yyaxis right
plot(x1, y1, 'Color','#0072BD' , 'Linewidth', 2)

hold on
set(gca, 'FontSize',15)
yyaxis left
title ('Annually binned Magnetization and Rainfall')
plot(x1, y2, 'Color',[0.85,0.33,0.10],  'Linewidth', 2)

%% plotting age model

figure
x1= age_mdl (:,2);
y1= age_mdl (:,1);
plot(x1, y1,'-o', 'Linewidth', 2);
set(gca, 'FontSize',15);
title ('Age Model');
ylabel ('Distance from top, mm');
xlabel ('Ur Th age');
%% O_18 regression with rain


%mdl=fitlm((binned_Bz(:,2)*1e8),binned_O_18(:,2))
mdl=fitlm(Oxygen_18_annual(:,2), (annual_precip(:,2)))
%mdl=fitlm(binned_O_18(:,2), (binned_Bz(:,2)*1e8))
figure
set(gca, 'FontSize',15);
title('Regression:O_18 and Total Magnetization')
xlabel('Bulk Magnetization')
ylabel('O-18')
plot(mdl)

%% C_13 regression with rain


mdl=fitlm(Carbon_13_annual(:,2), (annual_precip(:,2)))
figure
set(gca, 'FontSize',15);
title('Regression:O_18 and Total Magnetization')
xlabel('Bulk Magnetization')
ylabel('O-18')
plot(mdl)
%% Detrending carbon13 data

close all
figure
t1 = time_seriesCARBON;
A1 = x5;
D1 = detrend(A);

plot(t1,A1)
hold on
plot(t1,D1)
plot(t1,A1-D1,":k")
legend("Input carbon Data","Detrended Carbon 13 Data","c13 Trend")

%% Detrended O-18

figure
t2 = time_seriesOXYGEN;
A2 = x4;
D2 = detrend(A2);


plot(t2,A2)
hold on
plot(t2,D2)
plot(t2,A2-D2,":k")
legend("Input Oxygen Data","Detrended Oxygen Data","O18 Trend")

%% Detrend rain

close all
figure
t3 = time_seriesRAIN;
A3 = x3;
D3 = detrend(A3);

plot(t3,A3)
hold on
plot(t3,D3)
plot(t3,A3-D3,":k")
legend("Input Rain Data","Detrended Rain Data","Rain Trend")

%% Detrend limb mag
figure
t4 = (time_seriesALLLIMBS);
A4 = x1;
D4 = detrend(A4);


plot(t4,A4)
hold on
plot(t4,D4)
plot(t4,A4-D4,":k")
legend("Input Limb Bz Data","Detrended Limb Bz Data","Limb Bz Trend")

%% Regression centeral col to rain plots
%close all
%overlappyingarray matrix values
%column 1 is time; 
% column 2 is rain; 
% column 3 is limb Bz; 
% column 4 is center Bz.  All are smoothed.


mdl=fitlm((alldataoverlapping(:,2)/rainscale), alldataoverlapping(:,4))
figure
set(gca, 'FontSize',15);
plot(mdl)
title('Rain and Central Total Magnetization')
xlabel('central Magnetization')
ylabel('precip')

%% Regression centeral col to limb plots
%close all
%overlappyingarray matrix values
%column 1 is time; 
% column 2 is rain; 
% column 3 is limb Bz; 
% column 4 is center Bz.  All are smoothed.


mdl=fitlm(alldataoverlapping(:,3), alldataoverlapping(:,4))
figure
set(gca, 'FontSize',15);
title('Limb and Central Column Magnetization Regression')
xlabel('Limb Magnetization')
ylabel('Centeral Column Magnetization')
plot(mdl)

%% regression limb mag to rain

figure
mdl2=fitlm(alldataoverlapping(:,2), alldataoverlapping(:,3))
figure
set(gca, 'FontSize',15);
title('Rain and Total Magnetization')
xlabel('limb Magnetization')
ylabel('precip')
plot(mdl2)

%% Regression detrended data rain to o18


figure 
mdl3=fitlm((annual_precip(:,2)), D2)
figure
set(gca, 'FontSize',15);
title('Rain and Total Magnetization')
xlabel('limb Magnetization')
ylabel('precip')
plot(mdl3)

%% Regression rain to O18

figure 
mdl4=fitlm((annual_precip(:,2)), x4)
figure
set(gca, 'FontSize',15);
title('Rain and o18')
xlabel('o18')
ylabel('precip')
plot(mdl4)

%% checking fov 18 dates

figure
hold on
title('Magnetization')
xlabel('yrs')
ylabel('bz')
hold on
subplot(1, 2,1)
figure
y1=(Bzcropped_profile2(:,6));
x1= (Bzcropped_profile2(:,5));
plot(x1, y1, 'Color', '#D95319', 'Linewidth', 2)
hold on
subplot(1, 2,2)
x2= (Bzcropped_profile2(:,3));
set(gca, 'FontSize',15)
plot(x2, y1, 'Color', [0.2 .5 0.5], 'Linewidth', 2)

%% age model with error bars

figure

x = [1896 1963 2008 2016];
y = [127 53 10 1];
err = [2 2 2 4]
errorbar(x,y,err,'horizontal' , '-o','Color', [0.2 .5 0.5], 'Linewidth', 2)
set(gca, 'FontSize',15)

xlim([1894 2018])
ylim([0 131])
title('Age Model')
xlabel('Ur/Th Age')
ylabel('distance from top/mm')
legend

%% railsback

%plot O18 C rain and bz vertically for publication - 4 year
%close all
figure (1)


%1

y2= (time_seriesCENTER);
x2= (mag_interpCENTER);
set(gca, 'FontSize',15)
plot(x2, y2, 'Color', [0.2 .5 0.5], 'Linewidth', 2)


set(gca, 'FontSize',15)

ylim([1913 2016])
title('Magnetization')
xlabel('Bz')
ylabel('years')


%%
mdl9=fitlm(yearly_compiled_Bz(:,1), yearly_compiled_Bz(:,2))
figure
set(gca, 'FontSize',15);
title('cent vs flank')
xlabel('flank Magnetization')
ylabel('Centeral Column Magnetization')
plot(mdl9)

%% plots for Tyler, annual binned and raw - both unsmoothed
close all

figure
hold on
subplot(2,1,1)
x=(Bzcropped_profile2(:,5));
y1=abs(Bzcropped_profile2(:,2));



plot(x, y1,'Linewidth', 2,'Color', '#D95319')

xlabel('Calendar Year')
ylabel('bz raw')
xlim([1913 2016])


title('bz raw Center Column')

subplot(2,1,2); 
x3 = (yearly_compiled_Bz(:,1));
%y3 = (abs(yearly_compiled_Bz(:,2)*1e8));

y3 = (abs(yearly_compiled_Bz(:,2)));
yyaxis right
plot(x3,y3, 'Linewidth', 2)
title('Center column Bz 1 year bin')
xlim([1913 2016])
% ylim([0 1])

xlabel('Years')
ylabel('bz 1 year')
hold on
x4=rain_total(:,1);
y4= rain_total(:,2);
yyaxis left
ylabel('rain')
plot(x4,y4, 'Linewidth', 2)
xlim([1913 2016])
ylim([200 2300])




%% plot 1 year bin bz to rain unsmoothed

figure
x1=yearly_Bz_rain(:,1);
y1=(abs(yearly_Bz_rain (:,2)));
y2=yearly_Bz_rain (:,3);
yyaxis right
plot(x1,y1)
title ('Bz1 yr')
hold on
yyaxis left
plot(x1,y2)
xlim([1913 2016])

%%  regression 1 year bin bz to rain unsmoothed
mdl10=fitlm(yearly_Bz_rain (:,3), yearly_Bz_rain (:,2))
figure
set(gca, 'FontSize',15);
title('annual raw precip to annual cent bz not smoothed')
xlabel('rain')
ylabel('Centeral Column Magnetization')
plot(mdl10)

%%  figure plot of limb raw/annual/data for supplemental






all_limbs_raw={mag_interpAB_raw, mag_interpBC_raw, mag_interpCD_raw,mag_interpDE_raw};

% 
% if plot_smoothed
%     figure
%     plot(time_seriesBC, mag_interpBC_raw)
%     hold on
%     plot(time_seriesAB, mag_interpAB_raw,'r')
%     hold on
%     plot(time_seriesCD, mag_interpCD_raw,'g')
%     hold on
%     plot(time_seriesDE, mag_interpDE_raw,'m')
%     hold on
%     plot(time_seriesCENTER, -1*mag_interpCENTER,'k--')
%     hold on
%     plot(time_seriesRAIN, rainscale*rain_interpRAIN2,'b--')
%     hold on
%     plot(time_seriesOXYGEN, Oxygen_interp,'r--')
%     hold on 
%     plot(time_seriesCARBON, Carbon_interp,'g--')
% end
% 
% 
%%  spearman correlation coefficient for 1. rain vs center 2 rain vs limb and 3 center vs flank
%column 1 is time; 
% column 2 is rain; 
% column 3 is limb Bz; 
% column 4 is center Bz.  All are smoothed.

a = (alldataoverlapping(:,1)); %time
b = (alldataoverlapping(:,2)); %rain
c = (alldataoverlapping(:,3)); %flank Bz
d = (alldataoverlapping(:,4)); %center Bz



[RHO,PVAL] = corr(b, d,'Type','Spearman'); %rain & center

[RHO,PVAL] = corr(b, c,'Type','Spearman'); %rain & center

[RHO,PVAL] = corr(c, d,'Type','Spearman'); %rain & center


